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Abstract 

There are many aspects to consider when designing a Rosenbrock-Wanner- 
Wolfbrandt (ROW) method for the numerical integration of ordinary differen- 
tial equations (ODE’s) solving initial value problems (IVP’s). The process can 
be simplified by constructing ROW methods around good Runge-Kutta (RK) 
methods. The formulation of a new, simple, embedded, third-order, ROW 
method demonstrates this design approach. 

This paper outlines a recipe used by the author for developing ROW methods. 
Being an engineer, my approach differs somewhat from that of the mathematician 
(e.g., Hairer and Wanner [8]). This engineer breaks down the complex problem or 
ROW construction into two simpler problems; whereas, the mathematician typically 
retains all complexity and solves a single problem. The objective, in either case, is 
to be able to develop efficient and robust formulae for integrating systems of ODE’s 
that possess stiff regions in their solution paths. 

Let us consider the numerical integration of an IVP described by the non-autonomous 
system of ODE’s 


— = y = f(x,y(x)) with y(x 0 )=yo, (!) 

ax 

where y is a vector composed of m dependent variables { y 1 y 2 ■■■ y m } T with x 
being the single independent variable. 


1 Design 

Selecting a good RK method first simplifies the overall process of ROW design. Specif- 
ically, this author selects embedded methods to facilitate run-time error analysis, de- 
veloped after the manner of Merson [10], but more commonly referred to as Fehlberg 

'Computational Materials Laboratory, MS 105-1, NASA Lewis Research Center, 21000 Brook- 
park Road, Cleveland, OH 44135-3191. Email: afreedQlerc.nasa.gov. 


1 



[3, 4] methods. Such methods advance a solution via the formulae 

S — 1 > 

Yn+l = Yn + h + 0(h p+1 ) 

s — 1 

y n +i = Yn + 

i=0 ' 

The order of accuracy belonging to approximation y (the solution advanced) is p , 
while the order of the embedded approximation y (used in error assessment via y — y) 
is q, where q ^ p. The summation limit s denotes the number of stages in the 
integrator. Subscripts n and n+1 signify the n th and (n+l) th integration steps. 
Implied in Eqn. (2) is the relationship x n+1 = x n + h, where h is the size of the 
integration step. Parameters c, and c* are the weights of integration, typically lying 
within the interval [0, 1] . 

The Runge-Kutta derivatives k, are evaluated according to 


i-i 

kj = y n +h Aijkj ') , (3) 

3 = 0 

where vector components a, are the quadrature points of integration, and matrix 
components Aij are the explicit coupling coefficients. Our intuition and experiences 
have taught us that these quadratures ought to he within the interval [0, 1]; in other 
words, all kj derivatives should be evaluated locally within the interval [x n . x n+ \] of 
integration. 

One does not have complete freedom in selecting the parameters of Eqns. (2 & 3). 
They are constrained by a set of equations known as order conditions (e.g., see Hairer, 
Nprsett and Wanner [7, pp. 143-155]). Even so, there remains an aspect of design in 
RK construction, this because there are an insufficient number of order conditions to 
uniquely determine all the unknowns. The author chooses to apply this flexibility to 
specify good quadratures and weights, in the sense of classical integration theory. An 
additional flexibility arises in the next stage of development, which the author has 
used to minimize the error constants. 

This completes my first stage in ROW design. 

The second step of ROW construction is to fabricate a one- parameter family (in 
d, the Calahan [1] stability parameter) of Rosenbrock [13] integrators of the Wanner- 
Wolfbrandt [9, 16] 1 type, i.e., a ROW method, holding the RK parameters fixed. 
Embedded ROW methods are described by the same formulae for the dependent 
variables y n+ i and y n+1 as are used by explicit RK methods, viz., Eqn. (2). What 
dis tinguish ROW integrators from RK integrators axe their derivative- like functions; 


1 Wanner was a visiting professor at Chalmers University of Technology in Goteborg, Sweden in 
1975-1976 where he lectured on stiff methods. Wolfbrandt was a student there at that time. 
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specifically, 


l—l / 

ki = f (xn+dih, y n + h-Yl ^ykj) + h f bi 

j =0 V y 3=0 


RK (explicit) part 
with semi-implicit coupling 


3=0 

semi-implicit part 


Bu—d'di. 


( 4 ) 


( 5 ) 


The evolution of state residing in k,, i.e., function f, is the same here as it is in the RK 
derivatives of Eqn. (3). It is the presence of the non-autonomous gradient din/ dx with 
quadrature-like coefficients b u and the Jacobian di n /dy with semi-implicit coupling 
coefficients B lj , that provide the additional contributions of a ROW method. Implicit 
RK methods do not possess the semi-implicit part of Eqn. (4); rather, they extend 
the upper limit of the summation (appearing in the argument of f) to i for the semi- 
implicit case, or to s - 1 for the fully implicit case. 

Like RK methods, one does not have complete freedom in selecting the additional 
parameters of Eqn. (4), viz., the b t , B tJ and d. They too must satisfy a set of order 
conditions (e.g., see Hairer and Wanner [8, pp. 110-127]), some of which will have 
been satisfied a priori when adopting this design approach because the Oii Aij., Ci and 
Cj are considered known and satisfy the RK order conditions. It is useful to express 
the bi and B tj as functions of d at this stage of development. For methods with 
four or more stages, there are insufficient order conditions to uniquely determine all 
remaining unknowns. This extra flexibility, when it occurs, has been used by the 
author to design for a more optimal performance by reducing the error constants. 

This completes my second stage in ROW design. 

The final aspect of constructing a ROW method is the selection of d. This pa- 
rameter governs the stability characteristics of the method. Also affecting stability 
are the orders of the method’s integrators and the number of stages the method con- 
tains (e.g., see Hairer and Wanner [8, pp. 103-107]). Ideally, one would like both 
integrators in an embedded ROW method to be L-stable, although seldom will this 
be possible. As a minimum requirement, both integrators should be no less than A- 
stable. Some stability has been sacrificed by designers “since, roughly, smaller values 
of d give better error constants but increase instability, there is no optimal choice 
which optimizes simultaneously error constant and stability” [9]. 

One may want to iterate on steps two and three to seek an optimal compromise 
between performance and stability. The library of IVP’s cataloged by Enright and 
Pryce [2] is a useful proving ground for this purpose. This author encourages devel- 
opers to utilize stiff and non-stiff problem sets when assessing ROW performance. 
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2 Implementation 

Direct implementation of the semi-implicit equations in (4) is not very practical. For- 
tunately, a change in variable introduced by Wolfbrandt [16, pp. 106-107] removes 
the diagonal contribution from the right-hand side of Eqn. (4) resulting in the trans- 
formed formulae 


Yn+1 = Yn + h £ X i*. + o(w +1 ) 
i=0 

5~1 

y n +i = y n + h^2 + o(h q+1 ) 

i=0 


(6) 


and 

I - h -jp) • Ki = {(xn+aih^n+h^aijKj) +hbi ^ fijjKj . (7) 

^ j—o 3=0 

These equations contain the change in variables: (3 = ^1 — B -1 , q = A • B -1 , 
X T = c T - B _1 and x T = c T • B -1 . Because of the constraints B tl = d > 0 V i and 
Bij — 0 V i < j, matrix B has an inverse. 

There are several aspects worthy of comment that help transform a ROW method 
into an efficient and robust algorithm suitable for code implementation, which we have 
discussed more fully in [5]. In summary: use LU-factorization with full pivoting to 
decompose the matrix [^1 — h(df n /dy)]; employ a local interpolate^ made popular by 
Shampine [15], to permit dense output without requiring step-size reductions; use a 
stable step-size controller, like the proportional-integral (PI) controller of Gustafsson 
et al. [6], to monitor step size; estimate the initial step-size internally; construct an 
error estimate using some blend, e.g., 50/50, of relative and absolute errors, where 
the weights of absolute error are dynamically updated; and finally, utilize numerical 
evaluations for the non-autonomous gradient and Jacobian, thereby reducing the 
potential for programming errors, plus significantly simplifying the client module one 
must write to make use of such integrators. When combined with a ROW method, 
these various procedures will produce an efficient and robust algorithm for numerical 
integration. 


3 Example 

The first stage in my recipe is to choose or construct an embedded RK method 
with good quadrature and weights of integration. As a means for demonstration, 
consider the simple, 3 stage, third-order method given in Table 1, which has good 
quadrature, i.e., {0 \ 1} T , plus Simpson’s weights of integration, viz., {| | |} T . 
The coefficients listed in this table satisfy the order conditions of a 3(2) RK method — 
a third-order solution with an embedded second-order integrator for error assessment. 
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Table 1: A 3(2) Runge-Kutta method. 


a i 

Aij 

0 

0 

1 

I 0 

2 

2 u 

l 

-12 0 

Cj 

0 1 0 


1 2 1 

c j 

6 3 6 


Table 2: A 3(2) Rosenbrock-Wanner-Wolfbrandt method. 


&i 


bt 

Pij 

0 

0 

d 

0 

1 

2 

1/2 d 0 

0 

-l/d 0 

1 

l/d 2/d 0 

-d 

-2/d —4(2 — 3d)/d(l — 2d) 0 

Xj 

l/d l/d 0 

Xj 

7/6 d 2(3 — bd)/3d(l — 2d) 1/6 d 


For a 3(2) ROW method in 3 stages, where the a, and A xj are specified a priori, 
the number of remaining order conditions happens to be one less than the number of 
unknowns to be solved for. Consequently, the b x and B X j can be uniquely expressed 
as a one-parameter family of Calahan’s [1] stability parameter d. Their values, after 
transformation, are listed in Table 2. This completes my second stage of ROW design. 

The final stage of design is to select a stability parameter d. This example is 
a rare exception where both integrators in a ROW method can be L-stable. A 
second-order integrator in 3 stages (the embedded integrator) will be L-stable if 
d e [0.18042531,2.18560010] (see Hairer and Wanner [8, Table 6.4]). A third-order 
integrator in 3 stages will be L-stable if d = 0.43586652, which is the reciprocal of 
the second zero in the third-order Laguerre polynomial; therefore, selecting 

d = 0.43586652 

could, arguably, complete the design. 

Not yet addressed is the issue of the error constants. Requiring the integrator to 
be at least A-stable leads to the constraint that d € [1/3, 1.06857902] (see Hairer and 
Wanner [8, Table 6.3]). The embedded integrator is L-stable over this region. One now 
seeks an answer to the question: Is there an A-stable value for d whose selection would 
result in a significant performance gain over its L-stable value, thereby warranting a 
compromise in the integrator’s stability properties? 

In most developments of ROW methods there are an insufficient number of order 
conditions to uniquely define the b x and B l3 as functions of d alone. When this is 
the case, the designer should make attempts to minimize the error constants. Error 
constants are nothing more than the order conditions for that order which is one 
greater than the actual integrator. For example, in the 3(2) ROW method of Table 
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2, the embedded integrator fails to satisfy the third-order order conditions, which are 
two in number, resulting in the two error constants 

= -1/12 and ef = -(1 - 6d + 6d 2 )/6 . (8) 

Likewise, the approximation advanced as the solution arises from a third-order inte- 
grator that fails to satisfy the fourth-order order conditions, which are four in number, 
resulting in the four error constants 

e ( x 3) = 0 4 3) = — d 2 /6(l - 2d) 

ef = (1 + 8d)/24 ef = - (1 - 12 d + 36 d 2 - 24d 3 )/24 

These are the error constants for the integrator of Table 2 that one would want to 
minimize in an attempt to optimize performance. Figure 1 presents the norm for 
these error constants, calculated according to 

iM p, ll = (£(4 p> ) ) ■ (io) 

' i ' 

Viewing Fig. 1, it is immediately apparent that there is no advantage to considering 
values for d larger than the L-stable value. Choosing d — 2/5 leads to a 29% improve- 
ment in the norm for the error constants with only a minimal sacrifice in stability, 
while choosing d — 1/3 gives a 48% improvement but places stability at the boundary 
of the A-stable region. The compromise value of d = 2/5 seems, at first glance, to be 
a good design selection, but as we’ll see, this is premature. 

3.1 Numerical Exercise 

As a demonstration of capability, the performance of this new, 3-stage, 3(2) ROW 
method is compared against the performances of two, 4-stage, ROW methods of 
higher order — the 3(4) method Wolfbrandt [16, 11], MROS3, and the 4(3) method 
of Shampine [14], S-ROW, which is the Rosenbrock integrator promoted by Press et 
a 1. [12, pp. 738-742]. These constitute tough competitors. Also compared are the 
performances of the RK methods contained within these three ROW methods. 

The Brusselator from chemical kinetics [7, pp. 115-116] serves as a good example 
problem, i.e., 

Vo = Co 4- ylyi - (ci + 1 )y 0 

V\ = c x yo - ylyi 

which, by varying the values for constants Co and c\, varies its stiffness properties. 
This is illustrated in Table 3, wherein A denotes an eigenvalue from the Jacobian of 
Eqn. 11. Case i has a periodic solution; it is not stiff. Cases ii, in and iv successively 
increase in stiffness, and have monotonic solutions over the interval x 6 [0, 100], which 
is the interval of integration for all reported results, with initial condition yo(0) = 1.5 
and j/x (0) = 3.1. 
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d 


Figure 1: Error constant norms belonging to second- and third-order integrators in 
3(2) ROW method of Table 2. 


Table 3: Stiffness Ratios of the Brusselator. 


case 

Co 

Cl 

^max/^min 

i 

1 

5 

7 

ii 

1 

50 

2500 

Hi 

1 

500 

250000 

iv 

1 

5 000 

25000000 
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Table 4: Performance of Three RK Methods. 


case 

RK 3(2) 

RK in MROS3 

RK in S-ROW 


e = 0.01 

i 

363/588 

491/569 

333/1 090 

a 

1955/2 441 

2 723/3 256 

980/1 129 

in 

19941/24 927 

32 451/163073 

9875/12 339 

iv 

198 791/248493 

321204/1616498 

98522/123151 


e = 0.001 

i 

638/1563 

696/1 650 

763/4883 

ii 

1969/1984 

2 675/2 762 

1004/1016 

Hi 

19 955/24 942 

31921/161011 

9899/12 370 

iv 

198 803/248 497 

311520/1420202 

98544/123174 


e = 0.0001 

i 

1 188/3 613 

1589/10889 

1881/16996 

ii 

1992/2 007 

2 715/2 773 

1060/1071 

Hi 

19979/24 965 

28370/73977 

9963/13 587 

iv 

198 828/248520 

316 566/1545963 

98602/123 235 


The performance of the three, contained, RK methods is tabulated in Table 4. 
Reported are the number of steps required over the number of attempts made, each 
at three different accuracies for step control, i.e., 0.01, 0.001 and 0.0001. RK 3(2) out 
performs the other two RK methods for the non-stiff ODE’s of case z, which should 
not be a surprise since neither of the other two constructions gave consideration to 
the RK method within. When the ODE’s of Eqn. (11) become stiff, the RK method 
in S-ROW is superior, although one would not normally use an RK method for such 
stiff problems. Accuracy of solution, e, has little influence on the solution time in any 
of the stiff cases. Stability controls the solution time here. 

The performance of the three ROW methods is listed in Table 5. Without ex- 
ception, all three ROW methods were more efficient then their embedded RK coun- 
terparts for all four cases, as measured by the number of integration steps required. 
There are slight improvements in the non-stiff case where accuracy controls step size, 
and dramatic improvements in the three stiff cases where stability controls step size, 
as one would expect for both situations. At the tighter tolerance of e = 0.0001, the 
fourth-order integrators MROS3 and S-ROW begin to out perform my third-order 
method, also an expected result. Since the 3(2) integrator of Table 2 has one less 
stage in it, it can make up for slight performance losses (as measured by the number 
of steps taken) due to the smaller system required for LU decomposition and the for- 
ward/backward substitution process accompanying it. The very slight performance 
gain obtained in an effort to minimize the error constants, as reported in Table 5, 
does not warrant compromising stability in this case, answering the question posed 
in the previous section. 
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Table 5: Performance of Three ROW Methods. 


case 

d= 1/3 

ROW 3(2) 
d = 2/5 

L-stable 

MROS3 

S-ROW 


e = 0.01 

i 

289/344 

249/367 

315/403 

231/268 

287/347 

a 

25/25 

25/25 

25/25 

21/21 

24/24 

Hi 

28/28 

27/27 

27/27 

24/24 

29/29 

iv 

30/30 

30/30 

30/30 

26/26 

31/31 


rb 

II 

O 

O 
o 
>— 1 

i 

481/880 

516/1220 

544/999 

324/424 

441/554 

ii 

36/36 

37/37 

37/37 

27/27 

34/34 

in 

40/40 

40/40 

41/41 

31/31 

39/39 

iv 

42/42 

43/43 

43/43 

33/33 

42/42 


e = 0.0001 

i 

897/2576 

976/3059 

1021/3 159 

494/632 

714/1318 

ii 

57/57 

59/59 

60/60 

38/38 

49/49 

Hi 

62/62 

64/64 

64/64 

43/43 

56/56 

iv 

65/65 

68/68 

68/68 

46/46 

61/61 


4 Remarks 

The author finds it much simpler to follow the design approach outlined in this paper 
when constructing ROW methods than to use the more brute force approach of his 
predecessors. 

The 3(2) ROW method presented herein is L-stable and computationally efficient. 
It was designed for use in simulation and optimization codes, where robust perfor- 
mance matters and the accuracy of solution can oftentimes be relaxed, hence the 
lower order. It is an excellent integrator for these purposes. 

Because of the linear system of equations that must be solved, more CPU work 
is required per step from a ROW method than from a RK method of comparable 
size. RK methods are therefore still preferred for non-stiff ODE’s. For stiff ODE’s, 
however, ROW methods are clearly superior. If one is uncertain about the stiffness 
properties belonging to a system of ODE’s, then selecting a ROW method should 
not excessively increase the cost of computation whenever the problem turns out to 
be non-stiff, but whenever it ends up being stiff, it will likely save on cost, and this 
savings can be significant. 
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